program define estimate_IRF_ML, eclass
	
	syntax varname, [L(int 8) Alpha(real 0.05)] 
	
	tempvar res
	
	* Create matrices for storing results
	matrix `res' = J(`l', 7, .)
	matrix colnames `res' = lag IRF se lb ub df_resid K
	
	* Estimate MA(l) model
	quietly arima `varlist', ma(1(1)`l')
	if e(converged) != 1 {
		display as error "ML estimation did not converge"
		exit 111
	}
	
	* Calculate degrees of freedom
	local df_r = e(N) - e(k)
	
	* Collect all results
	forvalues ll = 1/`l' {
		local t_crit = invt(`df_r', `alpha' / 2)
		local lag = `ll' 
		matrix `res'[`ll', 1] = `ll'
		matrix `res'[`ll', 2] = _b[ARMA:L`lag'.ma]
		matrix `res'[`ll', 3] = _se[ARMA:L`lag'.ma]
		matrix `res'[`ll', 4] = _b[ARMA:L`lag'.ma] + ///
			`t_crit' * _se[ARMA:L`lag'.ma]
		matrix `res'[`ll', 5] = _b[ARMA:L`lag'.ma] - ///
			`t_crit' * _se[ARMA:L`lag'.ma]
		matrix `res'[`ll', 6] = `df_r'
		matrix `res'[`ll', 7] = .
	}
	
	ereturn matrix s `res'
	
end 
